Quantum phase transition in a far from equilibrium steady state of XY spin chain 
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Using quantization in the Fock space of operators we compute the non-equilibrium steady state in 
an open Heisenberg XY spin 1/2 chain of finite but large size coupled to Markovian baths at its ends. 
Numerical and theoretical evidence is given for a far from equilibrium quantum phase transition 
with spontaneous emergence of long-range order in spin-spin correlation functions, characterized by 
a transition from saturation to linear growth with the size of the entanglement entropy in operator 
space. 
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Non-perturbative physics of many-body open quantum 
systems far from equilibrium is largely an unexplored 
field. In one-dimensional locally interacting quantum 
systems equilibrium phase transitions - quantum phase 
transitions (QPT) - can occur at zero temperature only 
and are by now well understood pQ. QPT are typically 
characterized by vanishing of the Hamiltonian's spectral 
gap in the thermodynamic limit at the critical point, 
and (logarithmic) enhancement of the entanglement en- 
tropy and other measures of quantum correlations in the 
ground state [2] . Much less is known about the physics of 
QPT out of equilibrium, studies of which have been usu- 
ally limited to near equilibrium regimes or using involved 
and approximate analytical techniques (e.g. [3j[4]). 

There exist two general theoretical approaches to a 
description of non-equilibrium open quantum systems, 
namely the non-equilibrium Green's function method [5], 
and the quantum master equation [6[[7]. In this Letter we 
adpot the latter and present a quasi-exactly solvable ex- 
ample of an open Heisenberg XY spin 1/2 chain exhibit- 
ing a novel type of phase transition far from equilibrium; 
characterized by a sudden appearance of long-range mag- 
netic order in non- equilibrium steady state (NESS) as the 
magnetic field is reduced, and the transition from satu- 
ration to linear growth with size of the operator space 
entanglement entropy (OSEE) of NESS. 

The Hamiltonian of the quantum XY chain reads 



and the last spin are chosen (M = 4) 
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where cr%i y,z , m = 1, . . . , n are Pauli operators acting on 
a string of n spins. We may assume that parameters 7 
(anisotropy) and h (magnetic field) are non-negative. It is 
known that XY model ([I]) exhibits (equilibrium) critical 
behavior in the thermodynamic limit n — > 00 along the 
lines: 7 = 0, h < 1, and h = 1. Here we consider an open 
XY chain whose density matrix evolution p(t) is governed 
by the Lindblad master equation [6 (we set ft = 1) 



dp 

dt 



Cp := -i[H, p] + J2 ( 2L ^ L l ~ i L l L ^ P}) ( 2 ) 



and study a phase transition in NESS. The simplest non- 
trivial bath (Lindblad) operators acting only on the first 
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where = (cr^ ± icr^)/2 [8]. For h ^> 1, the ratios 
r^/r^ = exp(— 2h/T\) are simply related to canonical 
temperatures of the end spins T\, A = L, R. 

Note that Lindblad equation Q can be rigorously 
derived within the so-called Markov approximation [7 
which is justified for macroscopic baths with fast inter- 
nal relaxation times. As shown in [9 , Eq. ([2| with 
( Tp ) can be solved exactly in terms of normal master 
modes (NMM) which are obtained from diagonalization 
of An x An matrix A written in terms of 4 x 4 blocks 



A z , m = £/, m (-2/iR + £ U B L + <^ n B R ) 
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where R 7 

T\)(jy <g> 1 2 + ^ - r$)(a z + ia x ) ® a» 

Following [9], the key concept is 4 n dimensional Fock 
space of operators JC spanned by an orthonormal ba- 
sis Pry 



a-i ( a 2 « 2) , • — u i a '2 ' ' ' LL '2n ' lx j ^ {0, 1} where 

W 2 m-1 = 0-^U rn , <rn CT Z rn ,,W 2 rn = ^Um^m^m' are 2fl 

anticommuting Majorana operators {vjj,Wk} = 2Sj,k- 
We introduce canonical adjoint Fermi maps over /C, de- 
fined as Cj\P a ) = S aji i\wjP a ), so the quantum Liouvil- 
lean ^ becomes bilinear C = a • Aa + const 1 in Her- 
mitian maps a 2 j-i = ~^(cj + c]),a 2 j = ~^{cj - ct), 
satisfying {a pi a q } = S p ^ q . Note that the eigenvalues of 
An x An antisymmetric matrix A Q called rapidities come 
in pairs ft, -ft, ft, -ft, ••• ,ftn ) -ftn, Reft > 0. The 
corresponding eigenvectors v p ,p = 1, . . . , 4n, defined by 
A^2j-i = P3H23-11 ^^2j = —PjV.2ji can always be nor- 
malized as v_2j-i ' v 2 j = 1 and v p -y_ q = otherwise. 

Writing NMM maps as bj — v 2 j_ 1 -a, b f - = v 2 j -a, in gen- 
eral bj bp obeying almost-canonical anticommutations 
\bj,b k } = {bpb' k } = 0, {bj,b f k } = the Liouvillean 
(ij) takes the normal form, t = —^Y^jZiPj^jh- Thus 
a complete set of A n eigenvalues of C (real parts being 
the relaxation rates) can be constructed as — 2^- z^-ft 
where Vj G {0, 1} are eigenvalues of 2n mutually com- 
muting, non-hermitian number operators bjbj. 
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FIG. 1: Rapidity spectrum {(3j} around the imaginary axis 
for n = 640, = 0.5, = 0.3, Tf = 0.5, = 0.1, 7 = 0.5, 
and h = 0.3 < h c (left, blue) and h = 0.9 > h c (right, red), 
compared to dispersion |7| of the free XY model (center), 
where dashed curve indicates the critical case h — h c = 0.75. 

Let |ness) be the element of JC corresponding to the 
stationary solution p ness (NESS) of Eq. ([2|, i.e. zero 
eigenvalue of £, Vj = 0. The main result of [9] (Th. 3) 
takes into account the fact that |ness) is a right- vacuum 
of C - the left-vacuum being the trivial identity-state 
|1) - and asserts that any quadratic physical observable 
can be explicitly computed in terms of eigenvectors v , 
tr (wjW k p ness ) = S jyk + (l\cjC k \ness), 

^ 2n 

(l\CjC k \neSs) = - ^(^2m,2j-1^2m-l,2fc-l ~ (5) 
m=l 

— V2m,2jV2m-l,2k—iV2m,2jV2m-l,2k-l—iV2m,2j-lV2m-l,2$ • 

Higher order observables can be computed using the 
Wick theorem. For example, noting a z m = — i^2m-i^2m, 
spin-spin correlator which we shall study later reads 

Cz, m = tr (erf cr^p ness ) - tr (erf p ness ) tr (cr^p ness ) (6) 
= ( 1 1 C21- 1 c 2m I ness) (1 1 c 2 i c 2m - 1 1 ness) 
- (l|c2z-iC2 m -i|ness)(l|c2^c 2m |ness) if I ^ m. 

As proven in [9 j, NESS is unique iff rapidity spectrum is 
non-degenerate, (3j / for all j, and (almost) any initial 
state aproaches NESS asymptotically exponentially with 
the rate A = 2 mim,- Re if A > 0. 

Let us now proceed to detailed analytical and numer- 
ical investigation of the structure of NESS in XY chain. 
The bulk spectrum of rapidities for n — > 00 is insensitive 
to the coupling to the baths and is given by (3 = ±ie(0), 
(j) G (— 7r, 7r] where 



e((j>) = y (cos (j) — h) 2 + 7 2 sin 2 (j) (7) 

is the quasi-particle dispersion relation in an infinite XY 
chain (see e.g. [10]). For a finite chain ([T]) with the bath 
coupling on the edges ([3| we find that the bulk (nearly 
continuous) rapidity spectrum gains a small never van- 
ishing real part Re /3(<p) = 0(n _1 ). At the spectral edges 
/?*, P*\n=oo = ±ie(0*), with 0* defined by de(0*)/d^ = 
0, the gap is actually much smaller Re/3* = 0(n~ 3 ) (an- 
alytical result, generalizing [9]). Thus, the asymptotic 




FIG. 2: Spin-spin correlation matrices of NESS ^ for three 
different sizes n (columns) and different values of Ah — h — h c 
(rows) closely surrounding the critical value (pi. 7, same as 
in Fig.[l] Color-scale is proportional to log yOi,m\ and ranges 
from loglO -18 (dark-blue) to log 1 (white). 

relaxation time to NESS 1/A = 0(n 3 ) diverges in the 
thermodynamic limit n — > 00. 

We note, however, that the structure of the quasi- 
particle spectrum e(<j>) qualitatively changes as the mag- 
netic field crosses a critical value 

fc c (7) = 1 - 7 2 , (8) 

namely for h < h c the minimal quasi-particle energy 
exists for a nontrivial value of quasi- momentum (jf = 
arccos[ft/ft c (7)] yielding a new, non-trivial band edge (3*, 
whereas for h > h c the band edges can exist only at 
points (j)* = 0,7r (see Fig. [I]). Consequently, complex 
rapidities of an open XY chain shape up a third conden- 
sation point near the imaginary axis for h < h c which 
is composed of NMMs (eigenvectors of A) with pseudo- 
momenta near </>* 7^ 0, tt and has a dramatic effect on the 
structure of NESS as we demonstrate below. 

Indeed, as h < h c , we find the emergence of long 
range magnetic correlations (LRMC) characterized by 
non-decaying structures in the correlation matrix C^ m 
Typical size £ of the correlation patches is of the or- 
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FIG. 3: Spin-spin correlator C(r) computed as an average of 
Cz, m with fixed r = \l — m\ and \l + m — n\ < 0.08n. In (a) 
we plot C(r) in the critical case h — h c — 0.75 (7, same as 
in Fig. [I]) for several sizes n — 160, 320, 640 (indicated), while 
dashed line indicates asymptotic r~ 4 decay (double log-scale). 
Inset shows the scaled correlator n u C{r) versus r/n with v — 
4.09 for the same data (normal- log scale). In (b) we plot C(r) 
for changing h = 0.7505,0.751,0.752,0.755,0.76,0.77 > h c 
(right-to-left colored curves) in (log-normal scale) indicating 
exponential decay. Inset shows numerically determined local- 
ization length £ versus Ah — h — h c (points) as compared to 
theoretical estimate (|9| (dashed line). 

der £ ~ 1/0* (Fig.[2|. For h « h c one finds critical scaling 
0* « [2(h c — h)/h^/ 2 which agrees with the data. 

In the critical case h = h c (see Fig. [3^1) one finds power- 
law decay of the correlation matrix C^ m oc \l — m\~ 4 if 
neglecting finite size/boundary effects. If we scale the 
distance we find numerically a finite size scaling n u Ci jin = 
f(\l — m\/n) where v = 4.09 and f(x) is some function 
describing data for all large n (inset of Fig. [3^l) . Critical 
point h = h c is also characterized by faster closing of 
the spectral gap A of Liouvillean, namely there we find 
Re/3* = 0(n -5 ), meaning n 2 times longer relaxation 
times of generic solutions p(t) of 

For h > h c , we have 0* = and no LRMC in NESS. 
Then one finds an exponential decay of the correlation 
matrix C^ m oc exp(— \l — ra|/£) with the localization 
length which can be estimated theoretically from a scat- 
tering problem defined by the matrix Q: 

r 1 = ^ coslC 1 {h /h c ) « 4[2(/i - h c )/h c )]^ 2 (9) 

where factor 4 reflects the fact that C/ ?m is a 4-point 
function in NMM amplitudes v p (see Fig. [3]o) . 

The above results are summarized in a nonequilibrium 
phase diagram of XY chain (Fig. [4| showing the residual 

correlator C res = El!~ m|>n/2 <WEl!~ m|> " /2 1 (which 
is found to be always negative) in the j-h plane, with the 
critical curve h c (j) separating the two phases. Note that 
the other boundary lines 7 = (XX chain) and h = 
(XY with zero field) are not in LRMC phase. 

In analogy to equilibrium QPTs [TTJ [T2] . we wish 
to characterize the non-equilibrium transition in terms 
of quantum information theoretic concept, namely with 
the difficulty of classical simulation of p n ess which is de- 
scribed in terms of OSEE [13] (or block-entropy in /C), 




FIG. 4: Phase diagram showing the residual correlator — C res 
(log-scale indicated) calculated on 200 x 200 grid of parameter 
space h, 7 for size n — 160 (T^ same as in Fig. [lj. Red curve 
is the critical line h c = 1 — 7 2 . Note that C res is practically 
insensitive to increasing the size n in LRMC phase (bright). 

i.e. von Neumann entropy S(n) = — tr [ ljn / 2 ] Alog 2 R 
of the reduced density matrix of a half- chain R = 
tr [ n /2+i, n ] |ness)(ness|. tr corresponds to a partial 
trace over the sublattice [7, k] . Straightforward cal- 
culation, combining Refs. [9j [12], results in S(n) = 

-e;=i((| + %)iog 2 (i + + (l - %-)iog 2 (i - vj)), 

where rjj are n positive eigenvalues of an upper-left (or 
lower-right) 2n x 2n [14] block of 4n x An Hermitian corre- 
lation matrix ~D p , q = (ness|a p ag|ness)/(ness|ness). D can 
be computed by expressing a p in terms of NMM maps bj 
and b) (not ty), a = Q*6 + Qb\ Namely, D = Q*TQ T , 

where Tj^ = Yl^Li v< 2j-i,p v 2k-i p is 2nx 2n matrix, and 
Q = V G Ki2 where K12 designates upper-right 2n x 2n 
quarter of An x 4n matrix K = -(V | -V*)" 1 ^ -V*) 
and (V e ) p ,fc = v 2 k, P , (V )p )fe = v 2 k-i, P are 4n x 2n matri- 
ces. (X| Y) denotes vertical concatenation of two 4n x 2n 
matrices into a single 4n x 4n matrix. 

The resulting behaviour of S(n) in NESS of XY chain 
is striking (see Fig.j5|: LRMC phase h < h c is character- 
ized with a linear growth S(ri) = sn + const, with some 
constant s > 0. This has to be contrasted with a logn 
growth found for equilibrium critical models [12]. As h 
approaches h c the slope s approaches as s oc (h c — h) T ', 
with numerically determined critical exponent r w 0.80, 
and the fluctuations of 5(n) around an average linear 
growth increase. These fluctuations can be explained by 
sensitive dependence of NESS on boundary conditions 
(bath couplings or size changes) due to long range corre- 
lations, evident also in the structures of the correlation 
matrices (Fig. [2j. Note also an interesting 'quantiza- 
tion of bipartite entanglement' which is observed for very 
small h c — h where S(n) can take only approximately a 
discrete set of values S(n) ~ So + fc, k G Z + and which 
can be explained by the quasi-particle picture of NMM. 
At and above the critical field h > h c we find satura- 
tion S(n) = O(l), and vanishing fluctuations of S(n) 
since there NESS becomes insensitive to boundary con- 
ditions due to fast decay of magnetic correlations. Only 
there can NESS be efficiently simulated, e.g. in terms 
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FIG. 5: OSEE S(n) (symmetric chain bipartition) for 7, 
of Fig. [I] and different h < h c — 0.75 (indicated). Best fit- 
ting linear growths are indicated with straight lines. Dashed 
horizontal lines indicate So + 1 and 5*0 + 2, So being the sat- 
uration value for h = h c . Left inset just magnifies the scale 
while right inset shows the slope of S(n) growth vs. h c — h 
(log- log) and dashed line indicates \h c — h\ ' 8 . 

of matrix product states [15], by numerical methods like 
density matrix renormalization group (DMRG) [16]. 

All the numerical results presented above have been 
obtained for a fixed non-equilibrium bath couplings T\ = 
0.5, = 0.3, Tf = 0.5, = 0.1. However, the re- 
sults did not change qualitatively, in particular the phase 
boundary, when we (i) varied the bath couplings [8], 
(ii) coupled several spins around each end to Lindbla- 
dian baths, or (iii) even set the bath couplings equal 
= r^. The latter case (iii) does not represent 
equilibrium situation, i.e. p n ess is not a thermal state 
Pt = Z~ x exp(— H/T) as the XY chain is not ergodic 



[TP] , For example, no discontinuity at h = h c appears 
in the properties of pr for any T, and correlator C(r) 
essentially always decays with T— dependent rates [TP] , 
whereas in non-LRMC phase of NESS decay length £ is 
asymptotically insensitive to bath parameters Fur- 
thermore, thermal states in one-dimension have always 
bounded (in n) OSEE [17], and related quantities like 
mutual information [18], hence the simulation complex- 
ity of NESS is qualitatively different. 

In spite of demonstrated discontinuity in the spin-spin 
correlation function, the local observables such as en- 
ergy or spin density in NESS are numerically found to be 
smooth functions of h at h c , so the non-equilibrium tran- 
sition appears to be of high or infinite order (similar to 
Kosterlitz-Thouless transition). LRMC phase could per- 
haps be difficult to detect experimentally as the residual 
correlation C res is not larger than few times 10 -4 (Fig.[4| 
even in the optimal case (w.r.t. varying T^). 

In conclusion, we report on the QPT in NESS of open 
quantum XY spin chain, whose theoretical and numerical 
description is formally analogous to equilibrium QPTs in 
spin chains at zero temperature inasmuch as NESS can 
formally be treated as a 'ground state' of the quantum Li- 
ouvillean. We show that the phase transition is of mean- 
field type as the quasi-particle picture gives a satisfactory 
theoretical description, in particular the phase boundary 
between long-range and exponentially decaying magnetic 
correlations. We have demonstrated that the two phases, 
respectively, correspond to linearly growing and saturat- 
ing entanglement entropy of NESS in operator space as a 
function of the chain length. This behavior is drastically 
different than in equilibrium quantum XY chains. We 
thank M. Znidaric for useful comments and independent 
verifications of the results on small systems with DMRG 
[T6] codes. The work is supported by grants P 1-0044 and 
Jl-7347 of Slovenian Research Agency. 
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